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Abstract: In this work, we introduce GRChombo; a new numerical relativity code 
which incorporates full adaptive mesh refinement (AMR) using block structured Berger- 
Rigoutsos grid generation. The code supports non-trivial “many-boxes-in-many-boxes” 
mesh hierarchies and massive parallelism through the Message Passing Interface (MPI). 
GRChombo evolves the Einstein equation using the standard BSSN formalism, with an 
option to turn on CCZ4 constraint damping if required. The AMR capability permits 
the study of a range of new physics which has previously been computationally infeasible 
in a full 3 + 1 setting, whilst also significantly simplifying the process of setting up the 
mesh for these problems. We show that GRChombo can stably and accurately evolve 
standard spacetimes such as binary black hole mergers and scalar collapses into black 
holes, demonstrate the performance characteristics of our code, and discuss various 
physics problems which stand to benefit from the AMR technique. 
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1. Introduction 

Almost a hundred years after Einstein wrote down the equations of General Relativ¬ 
ity [1], solutions of the Einstein equation remain notoriously difficult to hnd beyond 
those which exhibit signihcant symmetries. Even for these highly symmetric solutions, 
basic questions remain unanswered. A famous example is the question of the non- 
perturbative stability of the Kerr solution - more than 50 years after its discovery, it is 
not known whether the exterior Kerr solution is stable. The main difficulty of solving 
the Einstein equation is its non-linearity, which dehes perturbative approaches. 

One of the main approaches in our hunt for solutions is the use of numerical 
methods. Numerical methods have been used to solve the Einstein equation for many 
decades, but the past decade has seen tremendous advances. A particular watershed 
moment was the breakthrough in evolving the inspiral mergers of two black holes [2-4] 
in 2005, a crucial milestone in the growth of numerical relativity as a discipline and as 
a tool. The other driver of this development is an explosion in the availability of large 
and powerful supercomputing clusters and the maturity of parallel processing technol¬ 
ogy such as the Message Passing Interface (MPI) and OpenMP, which open up new 
computational approaches to solving the Einstein equation. 

We anticipate that this development will continue to accelerate, partly driven by 
upcoming observational projects. The gravitational wave detector LIGO is expected to 
start Advanced LIGO science runs in late 2015, and there are hopes that the sensitivity 
might be good enough to achieve a hrst detection of gravitational waves from binaries. 
In the longer term, the European Space Agency (ESA) has designated the space-based 
eLISA detector an L3 launch slot (expected launch date around 2034), and the LISA 
Pathhnder spacecraft has a hrm launch date of late 2015. 

Beyond searching for gravitational waves and black holes, numerical relativity is 
now beginning to hnd uses in the investigation of other areas of fundamental physics. 
For example, standard GR codes are now being adapted to study modihed gravity [5], 
cosmology [6,7] and even string theory motivated scenarios [8-11]. In particular, there 
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is an increasing focus on solving GR coupled to matter equations in the strong-held 
regime: cosmic string evolution with GR, realistic black hole systems with accretion 
disks, non-perturbative systems in the early universe, etc. This nascent, but growing, 
interest in using numerical relativity as a mature scientihc tool to explore other broad 
areas of physics is one key motivation of this work. Since it is often difficult to have 
an intuitive picture of the entire evolution ahead of time, the code must be able to 
automatically adapt to ensure that all regions of interest always remain adequately 
resolved. 

In the numerical GR community, this requirement is largely met through a moving- 
box mesh rehnement scheme. This type of setup consists of hierarchies of boxes nested 
around some specihed centres, and the workhow typically requires the user to specify 
the exact size of these boxes beforehand. These boxes are then moved around, either 
along a prespecihed trajectory guided by prior estimates, or by automatically tracking 
certain quantities or features in the solution as it evolves. Boxes which come within a 
certain distances of each other may also be allowed to merge. A number of moving-box 
mesh rehnement codes have been made public over the recent years, many of which are 
built on top of the well-known CACTUS framework [12,13]. One such implementation is 
the McLachlan/Kranc code [14,15], which uses hnite difference discretisation and the 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) evolution scheme [16,17]. Similarly, the 
LEAN code [18,19], which uses the CACTUS framework, and BAM and AMSS-NGKU [20,21] 
also implement the BSSN formulation of the Einstein equations. There is also GRHydro 
which implements general-relativistic magnetohydrodynamics (MHD) for the Einstein 
Toolkit [22], building yet another layer of physics on top of evolution codes such as 
McLachlan/Kranc. There are also non-CACTUS codes such as SPeC [23] and bamps [24], 
which implement the generalised harmonic formulation of the Einstein equations using 
a pseudospectral method. In addition to these public codes, there is a plethora of 
closed-source codes. 

The moving-box mesh rehnement technique has found great success in astrophys- 
ically motivated problems such as two-body collision/inspiral. Outside of this realm, 
however, the setup can quickly become impractical, especially where one expects new 
length scales of interest to emerge dynamically over the course of the evolution. This 
can occur generically in highly nonlinear regimes, either by interaction between GR 
and various matter models, or by gravitational self-interaction itself which can exhibit 
complicated unstable behaviour in higher dimensions. In such situations, it is neces¬ 
sary to develop a code which has the hexibility to create rehnement regions of arbitrary 
shapes and sizes, anywhere in the computational domain as may be required. This can 
be achieved by using a fully adaptive mesh rehnement (AMR) technique, whose feature 
is generally characterised by the ability to monitor a chosen quantity at each time step 


- 3 - 



and insert higher resolution sub-regions where this quantity fails to he within some 
chosen bounds. Of course, the efficacy of such codes depend crucially on a sensible 
choice of these criteria, however when implemented correctly they can be an extremely 
powerful tool. The advantage here is twofold: AMR ensures that small emergent fea¬ 
tures remain well-resolved at all times, but also that only those regions which require 
this extra resolution gets rehned, thus allowing more problems to £t within a given 
memory footprint. To the best of our knowledge, PAMR/AMRD [25] and HAD [26] are the 
only two codes with full adaptive mesh rehnement (AMR) capabilities in numerical 
GR. 

In this work, we introduce GRChombo, a new multi-purpose numerical relativity 
code. GRChombo is built on top of the Chombo [27] framework. Chombo is a set of 
tools developed by Lawrence Berkeley National Laboratory for implementing block- 
structured AMR for solving partial differential equations. GRChombo features ^ include 
the following. 

• BSSN formalism with moving puncture-. GRChombo evolves the Einstein equation 
in the BSSN formalism. An option to turn on the CCZ4 constraint damping 
modification [28, 29] is also available. Singularities of black holes are managed 
using the moving puncture gauge conditions [3,4]. 

• Adaptive Mesh Refinement: Chombo provides full adaptive mesh refinement with 
non-trivial nesting topologies via the Berger-Rigoutsos block-structured adaptive 
mesh algorithm [30]. The user only needs to specify regridding criteria, and 
Chombo does the rest. Kreiss-Oliger dissipation is used to control errors, from 
both truncation and the interpolation associated with regridding. 

• MPI scalability: GRChombo inherits the parallel infrastructure of Chombo, with 
ability to scale efficiently to many thousands of CPU-cores per run. 

• Standardized Output and Visualization: GRChombo uses Chombo’s HDF5 output 
format, which is supported by many popular visualization tools such as Visit. 
In particular, the output hies can be used as input hies if one chooses to continue 
a previously stopped run - i.e. the output hies are also checkpoint hies. 

In this paper, we will detail these capabilities of GRChombo and illustrate how they 
expand the current held in numerical GR to permit new physics to be explored. The 

^Since the Chombo core is dimension-independent for up to six spatial dimensions, GRChombo could 
potentially be extended to simulate fully higher dimensional spacetimes without any symmetry as¬ 
sumptions, subject to computational resource availability. 
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design methodology, scaling properties and performance of GRChombo in a number of 
standard simulations are included. 

The paper is organized as follows: 

In Sec. 2 we describe the implementation of the code. In particular, in Sec. 2.1 
we establish the exact formulation of the equations which were used and our notation 
conventions, and in Sec. 2.2 we detail the AMR methodology and other “coding” 
aspects, such as hnite differencing, dissipation and load balancing. 

In Sec. 3 we give examples of several areas of physics which the code is well adapted 
to explore, and in which it offers advantages over existing codes. 

In Sec. 4, we present the results of standard tests, including the Apples with 
Apples tests [31], black holes and black hole mergers, and critical collapse. We test the 
AMR capabilities of the code, its robustness to regridding errors, and its scaling and 
convergence properties. 

We discuss our results and future directions in Sec. 5. 

Videos of several tests conducted in this paper, and examples of some new problems 
being tackled using the code, can be viewed via our website at http://grchombo.github.io 

2. GRChombo 

In this section, we will describe our numerical implementation of the Einstein equation 
in GRChombo. 

2.1 GRChombo equations and notation conventions 

The purpose of this subsection is to clearly state the equations of motion that we have 
implemented and our conventions for completeness. Since these are standard in the 
held, the experienced reader may want to skip this subsection. 

Many numerical relativity codes implement the so called BSSN formulation of the 
Einstein equation [16,17,32]. This formulation expresses the Einstein equation in a 
strongly hyperbolic form, and together with the “1 + log” slicing [33] and the “gamma- 
driver” gauge conditions [34], has allowed the stable simulation of dynamical spacetimes 
of interest, including black hole binaries. 

More recently, other rehned formulations of the Einstein equation based on the 
Z4 system [29,35] have been proposed, most notably the Z4c formulation [36] and the 
CCZ4 formulation [28].^ In the Z4 system, both the Hamiltonian and the momentum 

^Both the BSSN and CCZ4 equations have been written in a fully covariant form [37-39]. These 
covariant formulations can be advantageous in certain cases, and we plan to implement them in the 
future. 
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constraint are promoted to dynamical variables and hence constraint violating modes 
can propagate and eventually exit the computational domain. This may potentially 
result in a more stable evolution. In addition, the Z4 system can be augmented with 
damping terms so that constraint violating modes can be exponentially suppressed. In 
practical terms, the changes required between the CCZ4 equations and the standard 
BSSN equations are minimal and in GRChombo we have implemented both. 

In this work, we follow the indexing convention of [40]. The signature is (—h ++), 
and low-counting Latin indices a,b, ... are abstract tensor indices while Greek indices 
... denote spacetime component indices and run from 0,1, 2 , 3. Spatial component 
indices are labeled by high-counting Latin indices i, j, ... which runs from 1, 2, 3. Unless 
otherwise stated, we set G = 1 and c = 1. 

2.1.1 Evolution equations 

The Z4 system with constraint damping is [35] 

Rab + VaZb + VbZa-l^l [Ua Zb + Ub Za “ (1 + ^ 2 ) fi'afe^c] = 8 TT (Xab - ^ QabT^ (2.1) 

where Rab is the Ricci tensor associated with the metric g on the spacetime manifold 
Al, and V is the corresponding metric compatible covariant derivative. Tab is the 
stress-energy tensor of the matter and T = gab is its trace. If we set = 0, the Z4 
equations Eqn. (2.1) reduce to the standard (trace-reversed) Einstein equation. Here 
Ki and K 2 are parameters which control the damping. 

In the GRChombo code we use the standard 3-1-1 ADM decomposition of the space- 
time metric, 

ds^ = — dt^ + 'yij{dx^ + /3* dt)[dx^ -f (3^ dt ), (2.2) 

so that is the induced metric on the spatial slices and 

, (2.3) 

a 

is the corresponding timelike unit normal. The extrinsic curvature is dehned as 

Rij = —- {£n'y)ij , (2-4) 

where £ denotes the Lie derivative. As is customary, we decompose the induced metric 
as = p 7 jj so that det^j = 1 and y = (det 7 jj)~®. Similarly, the extrinsic curvature 
is decomposed into its trace, K = 7 *-^ Kij, and its traceless part so that 

Kij = — (^ij + 2 1 (2-5) 
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with 7 *-^ Aij = 0. In the Z4 system, one further dehnes 0 as the projection of the Z4 
four-vector along the normal timelike direction, 0 = —n^ Z^. Finally, the spacelike 
components of the four-vector, Zj, are included in a variable F* dehned as 

r = r + 27 *^'Zj, ( 2 . 6 ) 

where F* = and are the Christoffel symbols associated to the conformal 

metric 7 *^, 

^ + dk^ji - dnjk) . (2.7) 

Summarizing, the dynamical variables for the Z4 system are 

( 2 . 8 ) 

Setting to zero the Z4 four-vector, = 0, this system reduces to the standard BSSN 
system. 

Finally, we recall the various components of the matter stress tensor in the standard 
3-1-1 decomposition: 

p = nan,T-\ S, = -x^an,T-\ S,, = S = . (2.9) 

We are now ready to write down the evolution equations for CCZ4 system in the 
standard 3-1-1 decomposition [28]: 

dtX=^axK -^xdk/3’" + P^dkX, ( 2 . 10 ) 

dtlij = -2 a Aij -1- 7 ifc djP^ -f Xjk diP’" - ^ 7 ^ dkP’^ + dkpij , ( 2 . 11 ) 

dtK = -YWiDja + a{R + 2 DiZ^ + K^ -2X0)+ P^diK 

-3aKi{l +K2)e + 4:7ra{S-3p), (2.12) 

dtAij = yA \—DiDja -\- a {Rij + DiZj -\- DjZi — Stt a 5'jj)] 

+ (^Aij[K — 2 0 ) —2a Ail A^ + Aik djP^ + Ajk diP^ 

— 2 Aij dkP^ + P^ dkAij , (2.13) 

dte = ^a(^R + 2 D,Z" - Aij ^ - 2 0 

— Z'^ dia + P^ dkQ — aKi{2 + ^ 2)0 - Swap, (2-14) 

dtV = -2 dja + 2a A^^ - ^ f^djK - 3 

+ p^dkV + p^^d.dkP^ + ^ f^d.dkP'^ 
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+ ^ f* du(5’^ - + 2 «3 Q Z, dkl3^ - Zj dkl3^ 

+ 2 f djQ -Qdja-‘^aK Z^ - 2 a ki f - IGn a Sj . (2.15) 

Here Di is the metric compatible covariant derivative with respect to the physical metric 
'jij and [.. denotes the trace free part of the expression inside the parenthesis. The 
three-dimensional Ricci tensor, Rij, is split as 

Rij = Rij + I^j, (2.16) 


where 

+ f p,), -f 7'™(2ffpf,)fc^ + (2.17) 

1 _ . _ . 2 ~ ~ 

Rij = -{DiDjX + 'jijD^Dix) - -^ijD\Dix. (2.18) 

X X 

where Di is the metric compatible covariant derivative with respect to the conformal 
metric 7 p. Note that the three-dimensional Ricci Scalar is then R = x^'^Rij- 

Eqnations Eqn. (2.10)-Eqn. (2.15) are the CCZ4 evolntion eqnations as originally 
presented in [28], inclnding the extra damping parameter k^. This parameter controls 
the conpling of some qnadratic terms in the evolntion eqnation for T*. The choice 
K 3 = 1 corresponds to the fnlly covariant CCZ4 system, bnt as discussed in [28], 
it leads to instabilities in the evolution of spacetimes containing black holes. More 
recently, [41] showed that replacing ki —)■ Ki/a in Eqn. (2.10)-Eqn. (2.15) allows to 
stably evolve black hole spacetimes whilst retaining the full covariance of the CCZ4 
system. In GRChombo we have included a parameter that allows us to switch from the 
original formulation of the CCZ4 system to the more recent one proposed in [41], with 
the aforementioned redehnition of ki. 

Note that in the actual evolution, the values of the three-vector Zj are computed 
from the knowledge of the evolved variable T® and T®, which is computed from the 
conformal metric, Xij- Finally, we note that the evolution equations Eqn. (2.10)-Eqn. 
(2.15) reduce to the standard BSSN equations upon setting 0 = 0 and Z® = 0, and 
using the Hamiltonian constraint, Eqn. (2.23), in the evolution equation for K, Eqn. 
(2.12), to eliminate the Ricci scalar R. 


2.1.2 Gauge conditions 

To complete the set of evolution equations, we need to choose slicing conditions - we 
specify the gauge via driving conditions for the lapse a and shift /3® [34]. The optimal 
gauge conditions are in general physics dependent, and GRChombo allows the user to code 
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in whichever gauge conditions are well adapted to the application at hand. However, 
the most commonly used conditions which have been implemented in GRChombo are 
detailed below. 

The alpha-driver condition is usually written as a hrst order differential equation 

dta = + Hasfd^dia. (2.19) 

The commonly used 1 + log slicing applicable for black hole inspirals corresponds to 
fiai = 2, iia 2 = 1 and pag = 1. On the other hand, the maximal slicing condition, 
which preserves K = 0 and dtK = 0 at all slices, is a second order differential equation 

D'^a = a[KijK^^ + 47r(p + ^)], (2.20) 

which is useful for spherically symmetric collapse problems such as the critical scalar 
collapse scenarios. 

We specify the evolution equation for using the gamma-driver conditions [34], 

= 7]iB^ ( 2 . 21 ) 

dtB^ = - r]2B ^, ( 2 . 22 ) 

where is an auxiliary vector field, while pi, r] 2 -, and are input parameters. 
The usual hyperbolic gamma-driver condition uses the parameters rp = 3/4, = 1, 

/i /32 = 0 and 7^2 = 1. We have also included parameters that allow us to turn on 
standard advection terms in Eqn. (2.21)-Eqn. (2.22). 

In our tests in Sec. 4.2 and Sec. 4.4, where black holes are present, we manage the 
singularities with the so-called moving punctures method [3,4], which is a combination 
of the 1 -I- log slicing for a and gamma-driver for 13^. In addition, we hard code the 
condition a > 0 as is usual practice. 

2.1.3 Constraint eqnations 

GRChombo computes both the Hamiltonian constraint, 

H = R + K‘^ - KijK^^ - IQirp, (2.23) 

and the momentum constraint, 

M, = - diKji - - SttS,, (2.24) 

in order to monitor the accuracy of the calculation. 

For the algebraic constraints of BSSN, we do not enforce (by hand) the condition 
that the conformal metric has a determinant of one, but we do enforce after each 
timestep that Aij is traceless. 
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2.1.4 Scalar matter evolution equations 

We have included a single minimally coupled scalar field (j) as matter content 


with the equation of motion 




dV 

^ = 0 . 


d(j) 


(2.25) 


(2.26) 


As is usual, we decompose the second order Eqn. (2.26) into two first order variables 
0 and Hm 

= (2.27) 


a 


We note that our IIm is negative of fl in some references, e.g. [40]. Eqn. (2.26) is then 
decomposed into the following equations 


dtcf) = oRm + ld"di(p 


(2.28) 


and 


dtUM = milM + Y^{ad,d,cj) + d,<j)dia) + a (^KUm - . (2.29) 

We also use the energy momentum tensor of the scalar field 

V-^0 + 2V) (2.30) 

to calculate the matter components of the BSSN/CCZ4 system via Eqn. (2.9). 


2.2 GRChombo code implementation 

GRChombo is a physics engine built around the publicly-available adaptive-mesh frame¬ 
work Chombo [27]. GRChombo solves the system of hyperbolic partial differential equa¬ 
tions of the Einstein equation and scalar matter content (see section 2.1) using a finite 
difference scheme. 

A key feature of GRChombo is its highly flexible adaptive mesh refinement capability 
- to be precise, GRChombo uses Berger-Oliger style [42,43] adaptive mesh refinement 
with Berger-Rigoutsos [30] block-structured grid generation. GRChombo supports full 
non-trivial mesh topology - i.e. many-boxes-in-many-boxes. Morton ordering is used 
to map grid responsibility to neighbouring processors in order to optimize processor 
number scaling. 
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2.2.1 Discretization and Time-stepping 

We would like to evolve a set of fields in space (the state-vector <l>(x, t) = {0i, 02, 03, • • • }) 
through time t via the equations of motion 

( 9 $ 

(2.31) 

where is some operator on $ which, in the case of the Einstein equation, is non-linear. 

In GRChombo, both the space and time coordinates are discretized. Evolution in time 
is achieved through time-stepping t —)■ t -f At, where at each time step we compute 
the fluxes for each grid point individually. Time stepping is implemented using the 
standard 4th Order Runge-Kutta method, and hence, as usual, we only need to store 
the values of the state-vector at each time step. 

$ itself is discretized into a cell-centered grid. Spatial derivatives across grid points 
are computed using standard 4th order stencils for all spatial derivatives, except for 
advection terms which are implemented using an upwind stencil. The form of the 
stencils used exactly follows equations (2.2) through (2.6) of [44]. 

2.2.2 Berger-Rigoutsos Block-structured AMR 

GRChombo uses Chombo’s implementation of the Berger-Rigoutsos adaptive-mesh-rehnement 
algorithm [30], which is one of the standard block-structured AMR schemes. Block- 
structured AMR regrids by overlaying variable size boxes, instead of remeshing on a 
cell-by-cell basis (the “bottom-up” approach). The main challenge is to hnd an efficient 
algorithm to partition the cells which need regridding into rectangular “blocks”. In this 
section, we will briefly discuss the algorithm. 

Eor a given grid at some rehnement level I where / = 0 is the base level and Imax 
is some preset maximum rehnement level, we hrst “tag” cells for which rehning is 
required. The rehnement condition used by GRChombo is discussed later in this section. 
The primary problem of AMR is to efficiently partition this grid into regions which 
require adaptive remeshing. In block-structured AMR these regions are boxes in 3D 
or rectangles in 2D. Efficiency is measured by the ratio of tagged over untagged cell 
points in the hnal partitions. 

In each partition, we compute the signatures or traces of the tagging function 
f{x, y, z) of any given box 



1 f{x,y,z)dydz, 

(2.32) 

y(!/) = 

1 f{x,y,z)dxdz, 

(2.33) 
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f{x,y,z)dydx, 


(2.34) 



where f{x,y,z) = 1 if it is tagged for refinement and 0 otherwise. Given these traces, 
we can further compute the Laplacian of the traces d^X(x), dyY(y) and d^Z(z). Given 
the Laplacians, the algorithm can search for all (if any) inflection points individually 
for each direction - i.e. the locations of zero crossings of the Laplacian, and then pick 
the one whose 5{dfXi) is the greatest (corresponding to the line - or plane in 3D - 
separating the largest change in the Laplacian). This point then becomes the line of 
partition for this particular dimension. Roughly speaking, this line corresponds to an 
edge between tagged and untagged cells in the orthogonal directions of the signature. 
Furthermore, if there exists a point Xi with zero signature Xi{xi) = 0 (i.e. no cells 
tagged along the plane orthogonal to the direction), then this “hole” is chosen to be 
the line of partition instead. 

After a partitioning, we check whether or not each partition is efficient, specifically 
whether it passes a user-specihed threshold or fill factor, e < 1.0, 

Tagged Gells 
Total Gells 

If this is true, then we check if this box is properly nested^ [42,43] and if so we accept 
this partition and the partitioning for this particular box stops. If not, then we continue 
to partition this box recursively until either all boxes are accepted or partitioning no 
longer can be achieved (either by the lack of any tagged cells or reaching a preset limit 
on the number of partitions). Furthermore, GRChombo allows one to set the maximum 
partition size, which if exceeded will force a partitioning of the box. 

Note that a higher value of e means that the partitioning will be more aggressive 
which will lead to a higher efficiency in terms of hnal ratio of tagged to untagged 
cells - generating more boxes in the process. However, this is not necessarily always 
computationally better as partitioning requires computational overhead, which depends 
on the number and topology of the processors. The ideal £11 ratio is often a function 
of available processors, their topology and of course the physical problem in question. 

A partitioned box is then refined, i.e. its grids split into a finer mesh using the (user 
definable) refinement ratio = 5x^^^/5x\ and this process continues recursively until 
we either have no more tagged cells, or when we reached a preset number of refinement 

levels Imax- 



^Properly nested means that (1) a ^ + 1 level cell must be separated from an Z — 1 cell by at least 
a single I level cell and (2) the physical region corresponding to a / — 1 level cell must be completely 
filled by I cells if it is refined, or it is completely unrefined (i.e. there cannot be “half-refined” coarse 
cells). 
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Finally we need to specify a prescription for tagging which cells are required to be 
rehned. GRChombo tags a cell when any (set of) user selected helds 0 G $ pass a chosen 
threshold which sets a limit on the norm of the change in the value of the held 
across that cell, i.e. 


f{x,y,z) = < 


1 if \IY1 ^i) ~ ~ ^i))^ > ^(0) 

(2.36) 


0 


otherwise. 


This condition can be augmented, for example by using estimated truncation errors 
as tagging conditions instead. 

Partitioning can be done at every time-step for each rehnement level and this is 
a user preset choice per rehnement level. However, the user may wish to select a 
lower frequency because it might be useful to not partition at every timestep for a 
given rehnement level. One consideration is that it is important to let numerical errors 
dissipate (e.g. via Kreiss-Oliger dissipation, see Sec. 2.2.4) before remeshing. Once a 
new hierarchy of partitions is determined, we interpolate via linear interpolation from 
coarse to hne mesh, and average from hne to coarse mesh. 

Since the hner mesh has a smaller Courant number, each mesh level’s timestep is 
appropriately reduced via 

(2.37) 

GRChombo follows standard Berger-Collela AMR evolution algorithm [43]. Starting 
from the coarsest mesh, it advances the coarse mesh 1 time step i.e. t ^ t + Afh Then 
it advances the next hnest mesh n* times until the hne mesh “catches up” with the 
coarse mesh time. Once both coarse and hne mesh are at the same time t, GRChombo 
synchronizes them by averaging over the hne cells to the coarse cell values. We add 
that in a conservative system, this simple synchronization is not conservative and re¬ 
quires proper refluxing - the coarse huxes are replaced with a time-averaged hne mesh 
huxes. This step incurs additional overhead, and is at the moment not implemented by 
GRChombo as GR equations are not conservative. Nevertheless, we intend to implement 
conservative rehuxing as an option in a future version of GRChombo. 


2.2.3 Load Balancing 

GRChombo’s efficiency when running on a large number of distributed-memory nodes is 
highly dependent on efficient load balancing of the available computational work across 
those nodes. Load balancing seeks to avoid the situation where most of the nodes are 
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waiting for some small subset of nodes to finish their computational work, and it does 
this by seeking to distribute the amount of work to be done per time step evenly among 
all of the nodes. This can be non-trivial when AMR boxes at many different rehnement 
levels are simultaneously being evolved across the system. In addition, even within a 
single node, multiple OpenMP threads might be running, and the per-node workload 
needs to be balanced amongst those threads. 

For the inter-node load balancing, GRChombo leverages Chombo’s load balancing ca¬ 
pabilities to distribute the AMR boxes among the available nodes. It does this by 
building a graph of the boxes to be distributed, adding edges between neighbouring 
and overlapping boxes. A bin packing / knapsack algorithm is used to balance the 
computational work among nodes, where the work is assumed to be proportional to 
the number of grid points, and then an exchange phase is used to minimise the commu¬ 
nication cost. Because this load balancing procedure can be costly, we normally run it 
only every few time steps. In between runs of the load balancing procedure, new boxes 
generated by AMR rehnement stay on the node which holds the parent box. 

Within each node, the computational work is divided amongst the available OpenMP 
threads by iterating over the boxes to process using OpenMP’s dynamic scheduling ca¬ 
pability. This allows each thread to take the next available box from the queue of 
unprocessed boxes, instead of deciding ahead of time which boxes each thread will pro¬ 
cess. This is important because the boxes are varying in size. We generally divide even 
the coarsest level into multiple boxes so that it can be processed in parallel by multiple 
threads. 


2.2.4 Kreiss-Oliger Dissipation 

In a hnite difference scheme, instabilities can arise from the appearance of high fre¬ 
quency spurious modes. Furthermore, regridding generates errors an order higher 
than the typical error of the evolution operator, hence it is doubly crucial that we 
control these errors. The standard prescription to deal with this is to implement 
some form of numerical dissipation to damp out these modes. GRChombo implements 
N = 3 Kreiss-Oliger [45] dissipation. In this scheme, for all evolution variables u G 
{Aij, jij, iP, X, 0, F*}, the evolution equations are modified as follows 


dtUm t O^Uyti 


a 

64Aa; 


(^m-l-3 6 'U,jj_|_ 2 T 6'Um_2 Trim_3), 


(2.38) 


where min labels the grid point m, n the total offset from m and a is an adjustable 
dissipation parameter usually of the order This 3rd order scheme is accurate 

as long as the integration order of the finite difference scheme is 5 or less (which it is 
in our implementation using 4th order Runge-Kutta). 
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2.2.5 Boundary Conditions 

GRChombo supports both periodic (in any direction) boundary conditions, as well as 
any particular boundary conditions the user may want to specify (such as Neumann 
or Dirichlet types). A particular popular type of boundary condition is the so-called 
Sommerheld [34] boundary condition, where out-going radiation is dissipated away. For 
any held /, we impose the condition at the boundary 

df VXi df f - fo 

where r = ^Jx\ + X 2 + is the radial distance from the center of the grid, /o is the 
desired space-time at the boundary (typically Minkowski space for asymptotically hat 
spacetimes) and v the velocity of the “radiation”, which is typically chosen to be 1. 

2.2.6 Initial conditions 

GRChombo supports several ways of entering initial conditions. 

• Direct equations - Initial conditions which are described by known analytic equa¬ 
tions, such as the Schwarzchild solution, can be entered directly in equations form. 

• Checkpointing - The HDF5 format output hies from GRChombo doubles as check¬ 
pointing hies. A run can simply be continued from any previous state as long as 
its HDF5 output hie is available. 

• Entering from data hie - GRChombo allows one to insert data from a hie. 

• Relaxation - GRChombo has a rudimentary capability to solve for the initial metric 
given some initial mass distribution, and assuming a moment of time symmetry 
and conformal hatness. Given a guess for y, GRChombo relaxes it to the correct 
initial metric using a dissipation term which is proportional to a user chosen 
dissipation coefficient times the Hamiltonian constraint. 

The initial conditions used in this paper are mostly analytic or approximate an¬ 
alytic solutions, and so are entered directly into the code. In the critical collapse, a 
Mathematica numerical solution as a function of the radius is interpolated onto the 
initial grid. 

3. Using GRChombo for new physics 

The primary advantage of GRChombo over existing publicly available codes is its robust 
AMR ability. In this section we discuss several physical systems which can be studied 
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with GRChombo but that would be hard (if not impossible) to simulate using codes based 
on moving-box mesh rehnement. Most of the examples discussed in this section are 
still “work in progress” by the authors and the results will be presented in separate 
publications. 

3.1 Asymmetric scalar field bubbles 

One of the most fascinating and as yet not fully understood aspects of general relativity 
is the appearance of critical phenomenon in gravitational collapse as first discovered by 
Choptuik [46]. A comprehensive review can be found in [47]. 

Briefly, if we have an initial conhguration, such as a Gaussian shaped bubble of 
scalar held, and allow this to evolve under the action of gravity, the result will be either 
the formation of a black hole, or dispersal of the held to inhnity depending on the 
“strength” of the initial data. Varying any one initial parameter p of the conhgnration 
(snch as the height of the bubble), one hnds that there is a critical point p* at which 
the transition between the two end states occurs, and that the mass of the black hole 
created on the supercritical side follows the scaling relation, 

Moc{p-p*)\ (3.1) 

where the scaling constant 7 is nniversal in the sense that it does not depend on the 
choice of family of initial data. For a massless scalar in a spherically symmetric collapse, 
7 has been nnmerically determined to be around 0.37. 

The other key phenomenon which is observed is that of self-similarity in the solu¬ 
tions, or “scale-echoing”. Close to the critical point, and in the strong held region, the 
value of any gauge independent held 0 at a point x and time T exhibits the following 
scaling relation. 


(j){x,t) = (p{e‘^x, e^T), (3.2) 

where A is a dimensionless constant with another nnmerically determined valne of 
3.44 for a massless scalar held in the spherical case. The time T here is measured 
“backwards” - it is the diherence between the critical time at which the formation of 
the black hole occurs and the cnrrent time, with time being the proper time measured 
by a central observer. 

What one sees is therefore that, as the time nears the critical time by a factor of 
e‘^, the same held prohle is seen bnt on a scale smaller. This scale-echoing may be 
either continuous or discrete. In Choptuik’s seminal paper [46], a 1-1-1 adaptive mesh 
code was used to study such behaviour near the critical point. Since then there has 
been some progress in stndying the phenomenon in non spherically symmetric cases. 
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Figure 1: An asymmetric scalar field bubble collapse: In the left image, a slice through 
the centre of the field bubble is shown at a particular time during the collapse, with the 
instantaneous mesh superimposed. This same slice is shown in 3D in the right image to 
better illustrate how the mesh is adapted to the curvature. The AMR efficiently tracks the 
scales in the profile. In addition, new levels are added as the profile shrinks, leading to a 
more efficient computation in terms of memory and computational resources, and requiring 
less “human input” ahead of and during the run. 

see [24,48-52], but progress in making firm conclusions has been slower than expected, 
due to the extremely high rehnements required to study the stages of the collapse, 
which are magnihed three-fold in full 3-1-1 codes. 

We are currently investigating the problem of asymmetric bubble collapse with 
GRChombo. A snapshot of the evolution of one such example is shown in Figure 1. 
The mesh rehnes so that the held prohle is consistently resolved as the critical time 
is approached. Note that since the prohle is highly irregular, with a range of length 
scales represented in diherent (disconnected) regions, GRChombo provides a signihcant 
advantage over moving-box mesh rehnement in terms of computational cost. In addi¬ 
tion to the adaptation of the mesh to the local curvature, GRChombo can automatically 
add in new levels of rehnement during the evolution, which allows the scale-echoing be¬ 
haviour to be probed consistently as the prohle shrinks even in the absence of spherical 
symmetry. The results will be discussed in detail in a separate publication. 

3.2 Higher dimensional black holes/anti-de Sitter 

In recent years it has been realized that the dynamics of general relativity, even in 
vacuum, beyond the traditional asymptotically hat four-dimensional setting is much 
richer than previously anticipated. Some of these new directions involve considering 
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more than four spacetime dimensions and/or new boundary conditions, such anti-de 
Sitter (AdS) or Kaluza-Klein (KK) asymptotics. Black holes, as primary objects in any 
theory of gravity, in these new set ups exhibit two important new physical phenomena: 
Firstly, black holes with topologically non-spherical horizons are possible. The black 
ring of [53] is the hrst example of an asymptotically flat vacuum black hole with a non- 
spherical horizon. By now it is clear that this solution is the tip of the iceberg, and many 
more new types of black holes are known to exist. Secondly, vacuum black holes can 
be dynamically unstable under gravitational perturbations, as Gregory and Laflamme 
hrst showed in the case of black strings in asymptotically KK spacetimes [54] . Rapidly 
spinning asymptotically hat vacuum (and AdS) black holes suher from these Gregory- 
Lahlamme-type-of instabilities [55-57], and new types of non-axisymmetric instabilities 

[58] . In fact, anti-de Sitter itself is non-linearly unstable to the formation of black holes 

[59] , and the process is turbulent. There is a lot of interest in studying the dynamics of 
gravity, and black holes, in AdS motivated by the gauge/gravity correspondence [60]. 
See [61,62] for some (relatively recent) reviews with references. 

In a remarkable paper, [63] studied the endpoint of the Gregory-Lahamme insta¬ 
bility of black strings in hve-dimensions. This paper gave convincing evidence that 
the black string would pinch oh in hnite asymptotic time, thus providing a poten¬ 
tial counter-example of the weak cosmic censorship conjecture in non-asymptotically 
hat vacuum spacetimes. The evolution the Gregory-Lahamme instability for black 
strings showed that the horizon develops a fractal structure, with thin necks con¬ 
necting bulges at diherent scales. Moreover, the non-linear instability of AdS is of 
a turbulent nature [59], and in fact AdS black holes also suher from turbulent type- 
of-instabilities [64-67]. A common feature in all these instabilities is that in the fully 
non-linear regime, new length scales are dynamically generated and a priori one does 
not know where they will appear. Therefore, if one wants to use numerical GR to 
determine the endpoints of these instabilities, one needs a numerical method that can 
automatically resolve these newly generated length scales. Therefore, it seems that full 
AMR is not an option but a necessity."^ In addition, with GRChombo one should be able 
to simulate higher dimensional spacetimes with no symmetry assumptions using the 
fact that the Chombo core is dimension independent up to six spatial dimensions. This 
should hnd applications to the AdS/GFT correspondence, where one may be interested 
in simulating (4 -|- l)-dimensional asymptotically AdS spacetimes in full generality. 

Some of the aforementioned instabilities are currently being investigated using 
GRChombo and will be discussed elsewhere. In Fig. 2 we display snapshots of the 
meshes that are dynamically generated by GRChombo during evolution for the case of 

^The AMR capabilities of the PAMR/AMRD code were essential in [63]. 
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Figure 2: Snapshots of the rotation plane of dynamically generated meshes by GRChombo 
during the evolution of the instabilities of a black ring (left) and of a rapidly spinning spherical 
black hole (right). AMR ensures that the mesh is adapted to the non-trivial topology of the 
horizon and new levels are added where new structure appears. This is essential in order 
to have enough resolution where it is needed while keeping the computational cost of the 
simulation under control. 


an unstable black ring (left) and a higher dimensional black hole (right). GRChombo not 
only can adapt the mesh to the non-trivial topology of the horizon but it also adds new 
levels in regions where new structures appear during the evolution. This essential ca¬ 
pability ensures that all relevant length scales are correctly resolved whilst keeping the 
computational cost of the simulation under control. For these simulations, moving-box 
mesh rehnement would be prohibitively expensive and hence it is not a realistic option. 
Previous works on higher dimensional black hole physics in asymptotically flat spaces 
using moving boxes include [19,68]. In AdS, the code of [69] is based on PAMR/AMRD, 
whilst [9] uses a hxed domain decomposition and pseudospectral discretisation. 

3.3 A^-body problems 

Simulations of single black holes and binary mergers often make use of symmetries in 
the problem or Newtonian approximations to predict the levels of resolution required 
at each point in space, for each time-step. In three (or more) body problems, the 
trajectories of the objects are generally not known ahead of time and must be calculated 
numerically. The shift vector can be used to predict the movement of black hole centres 
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Figure 3: Triple black hole merger: Three black holes are evolved with GRChombo. The mesh 
is shown, which has adapted to the local curvature in Xj the variable plotted. 


locally, but this tracking must be added in for each black hole and the boxes adjusted 
accordingly. For GRChombo, it is trivial to add multiple black holes to a spacetime, 
without actively tracking their central points, and so these many-body systems are 
as easy to set up as binary ones (although they clearly require greater computational 
resource to run). For example, the GRChombo mesh can be seen adapting at each 
timestep in the triple black hole merger shown in Figure 3. 

Another (albeit rather contrived) example is the ring of black holes shown in Figure 
4. The set up of this ring was no more difficult than that of the binary or triple black 
hole spacetime - GRChombo automatically remeshes the grid given a set of analytic 
initial conditions without further user intervention, and there is no need to try and 
predict the paths of the black holes individually. It is clear that other fields of research, 
such as magnetohydrodynamics, would benefit greatly from this ability to maintain a 
consistent level of resolution throughout a simulation, and make this resolution follow 
the inherently unpredictable movement in a many body system. 

4. Testing GRChombo 

We detail the results of the standard Apples with Apples tests [31] in Sec. 4.1 when 
turning off AMR and using fixed resolution grids. In Sec. 4.2 we turn on the AMR 
abilities of the code and demonstrate that it can stably evolve spacetimes containing 
black-hole-type singularities. In 4.3 we study convergence of our code in a head on 
collision of two black holes. In Sec. 4.4 we demonstrate the ability of the code to evolve 
matter content by considering scalar helds with gravity, by recreating the results of the 
sub-critical and critical cases of Choptuik scalar field collapse detailed in [70]. In Sec. 
4.5 we discuss the weak scaling properties of the code on the Mira supercomputer. 
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Figure 4: Multiple Black Hole merger: A ring of black holes are evolved with GRChombo. 
The mesh is shown, which has adapted to the local curvature in x, the variable plotted. 


The test figures referred to in this section can be found in Appendix A. 

4.1 Apples with Apples tests 

In this section we describe the results of applying the code to the standard Apples with 
Apples tests in [31]. Here we give a brief description of the key features of the tests, 
but the reader should refer to the paper for full specifications. Where we do not specify 
details, our treatment can be assumed to follow that of the standard tests. The AMR 
capabilities of the code are not utilised in these tests, which were designed for a uniform 
resolution, in order to make our results comparable to other codes. (We consider the 
effects of regridding on code performance in Section 4.2.) 

4.1.1 Robust stability test 

The robust stability test introduces small amounts of random noise, scaled with the grid 
spacing, to all of the evolution variables, in order to test the code’s robustness against 
numerical noise. The test was conducted at resolutions of p = 4, p = 2 and p = 1, 
corresponding to grid spacings of 0.005, 0.01 and 0.02 respectively. No dissipation was 
added in the test. 

As shown in Figure 13, the error growth in the evolution variables did not increase 
with increasing grid resolution, and the Hamiltonian constraint H did not grow more 
for higher resolutions. Therefore, we conclude that the test is passed. 

4.1.2 Linear wave test 

A wave of fixed amplitude is propagated across the grid in the x-direction with periodic 
boundary conditions. The amplitude is small enough that the non-linear terms are 
below numerical precision, such that the behaviour under the Einstein equation is 
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approximately linear. The test measures the errors in magnitude and phase introduced 
by the code after 1000 crossing times. 

As can be seen from Figure 14, this error is 12 orders of magnitude smaller than 
the signal and therefore negligible. 

4.1.3 Gauge wave tests 

The BSSN formulation is known to produce unsatisfactory results for the gauge wave 
tests. GRChombo is no different in this respect. As can be seen in Figure 15, it becomes 
unstable after around 50 crossing times, with the Hamiltonian constraint increasing 
exponentially, even for a relatively small initial amplitude of the gauge wave of A = 0.1. 

As was shown in [28] stability can be achieved by adding in the CCZ4 constraint 
damping terms. GRChombo shows exactly this behaviour (hgure 15). 

4.1.4 Gowdy wave test 

The Gowdy wave evolves a strongly curved spacetime: an expanding vacuum universe 
containing a plane polarised gravitational wave propagating around a 3-torus. In the 
expanding direction we use the analytic gauge, dta = —The collapsing direction 
is evolved starting at f = to with harmonic slicing for the lapse and zero shift. A Kreiss- 
Oliger dissipation coeffecient of a = 0.05 was used in both directions. 

The results for both the BSSN and CCZ4 codes in the collapsing direction are 
shown in Figure 16, and in the expanding direction in Figure 17. 

As is found in the Apples with Apples tests [31] for other simple BSSN codes, 
GRChombo with BSSN and CCZ4 gives a less than satisfactory performance in this test in 
the expanding direction. The evolution is stable for approximately the hrst 30 crossing 
times, after which high frequency instabilities develop and cause code crash, due to the 
exponentially growing 'y^z component. In [31] it was found that this behaviour of BSSN 
could be controlled with dissipation, but that long term accuracy was not achievable. 

In the contracting direction the evolution is stable for the full 1000 crossing times 
and we were able to conhrm the convergence of our code. As shown in Figure 18, both 
BSSN and CCZ4 exhibit 4*^ order convergence initially. While convergence is never 
lost, the order is reduced at later times. This is similar to the behaviour found in [31] 
and [71]. 

4.2 Vacuum black hole spacetimes 

In this subsection we show that our code can stably evolve spacetimes containing black 
holes. 

All the simulations presented here used the BSSN formulation of the Einstein equa¬ 
tions, along with the gamma-driver and alpha-driver gauge conditions. Adding CCZ4 
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constraint damping gives better performance for the Hamiltonian constraint, as would 
be expected, but the results are broadly similar and so are not presented here. Unless 
otherwise stated, we perform the simulations with up to 8 levels of rehnement and we 
based our tagging/regridding criterion, Eqn. (2.36), on the value of y. We emphasise 
that the purpose of this subsection is to demonstrate that we can stably evolve black 
hole spacetimes, but we are not interested in extracting gravitational wave data or in 
studying convergence; this will be done in the next subsection. 

Where we refer to taking an norm of the Hamiltonian constraint H in a test, 
this is calculated as follows (using the weighted variable sum function in Visit): 



(4.1) 


where is Vi/Vtot, the fraction of the total grid volume V occupied by the ith box. 
Where the grid contains a black hole, we excise the interior by setting H to zero within 
the region in which the lapse a is less than 0.3 (which is an approximate rule of thumb 
for the location of an event horizon for a black hole in the moving puncture gauge). The 
difference in the results is small, since the error norm is dominated by regridding errors 
at the boundaries between meshes. We also exclude the values on the outer boundaries 
of the grid, which can distort the results in cases where periodic boundaries are used. 

4.2.1 Schwarzschild black hole 

First we evolve a standard Schwarzschild black hole in isotropic gauge, with a confor¬ 
mally flat metric, the lapse initially set to one everywhere, and the conformal factor 



(4.2) 


In this simulation, we chose the outer boundary of the domain to be at 600M and 
the spatial resolution in the coarsest mesh is lOM. We impose Sommerfeld boundary 
conditions. The initial value of x through a slice is shown in Figure 5. We see the 


expected “collapse of the lapse” at the singularity and the solution quickly stabilises 


into the “trumpet” puncture solution described in [72]. We hnd an apparent horizon 
and are able to evolve the black hole stably and without code crash for well over 
t = lOOOOM time steps as shown in Figure 19 (left). In this hgure we show the 
norm of the Hamiltonian constraint across the whole grid, and it remains bounded 
throughout the evolution. 

We monitor the ADM mass of the black hole by integrating over a surface near 
the asymptotically flat boundary, as seen in Figure 19 (right). We also monitor the 
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Figure 5: The profile for x through a slice perpendicular to the z axis is shown. 


angular momentum and linear momentum of the black hole, and find that these remain 
zero as expected, as shown in Figure 19 (right). These simple ADM measures rely on 
asymptotic flatness at the surface over which they are integrated, and so are sensitive to 
errors introduced by reflections at the boundaries, initial transients from approximate 
gauge choice or if the black hole is moving nearer the boundary (as in the boosted 
case). They are therefore less reliable as the simulation progresses, and we use them 
simply to conhrm that we are evolving the correct spacetime initially. 

4.2.2 Kerr black hole 

In this sub-subsection we present the results of a simulation of the Kerr black hole 
spacetime in quasi isotropic gauge as in [73] with the angular momentum parameter 
a = J/M set to 0.2. The domain size was chosen to be (320M)^ and the grid spacing 
in the coarsest level was AM. We impose periodic boundary conditions for simplicity, 
which limits the duration of the simulation due to boundary effects. 

In Figure 20 (left) we show the norm of the Hamiltonian constraint throughout 
the evolution. This plot shows that the amount of constraint violation remains stable 
during the simulation. In the right panel of Figure 20 we display the ADM measures 
for the three components of the angular momenta and the mass. This Figure shows 
that these quantities remain (approximately) constant during the simulation. 
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(b) Position at t = 100 



(c) Grid at t = 100 


Figure 6: Boosted black hole, movement: The boosted black hole moves across the grid 
diagonally with initial momenta of Px = 0.02, Py = 0.02 and Pz = 0.0, as expected, and the 
grid adapts to this movement, with the high resolution grids following the movement. 


4.2.3 Boosted black hole 

In this sub-subsection we evolve a boosted black hole using the perturbative approxi¬ 
mation from [40], with initial momenta set to Px = 0.02, Py = 0.02 and Pz = 0.0. The 
domain size was chosen to be (640)^, with spatial resolution in the coarsest grid of 4M. 
We imposed periodic boundary conditions at the outer boundaries of the domain. The 
black hole moves across the grid diagonally as expected, as is seen in Figure 6. 

In the left panel of Figure 21 we show the L? norm of the Hamiltonian constraint 
across the domain as a function of time. This plot shows that the constraints remain 
bounded throughout the simulation. In the right panel of Figure 21 we display the 
components of the ADM linear momentum during the simulation. In the continuum 
limit they should be constant and in our simulation they are indeed approximately 
constant. 

4.2.4 Binary inspiral 

In this sub-subsection we superpose the initial perturbative solution for two boosted 
black holes in [40], sufficiently separated, to simulate a binary inspiral merger. The 
domain size was (200M)^ with a grid spacing in the coarsest level of 5M. As in some 
of the previous tests, for simplicity we imposed periodic boundary conditions at the 
outer boundaries of the domain. 

We are able to evolve the merger stably such that the two black holes merge to 
form one with a mass approximately equal to the sum of the two. The progression of 
the merger is shown in Figure 7. The time evolution of the norm of the Hamiltonian 
constraint across the grid is shown in Figure 22. Again this remains stable throughout 
the simulation. 
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Figure 7: Binary Black Hole merger: Two black holes are evolved with GRChombo. The final 
stages of the merger are shown. 


4.3 Convergence test: head on collision of two black holes 

In this subsection we simulate the head on collision of two black holes and analyse 
the convergence of the code. We set up Brill-Lindquist initial data consisting of two 
static black holes of mass 0.5M with a separation of lOM, located at the centre of the 
computational domain. We extract the gravitational wave signal (see below). An initial 
burst of radiation is seen, which is a property of the superimposed initial data, prior to 
the main signal. Even though this set up could be simulated in axisymmetry, we have 
evolved the system without imposing any symmetry assumptions. So the results below 
correspond to a full 3 + 1 simulation with GRChombo. 

We performed runs at three different resolutions with 7 levels of refinement, each 
level having half the grid spacing as the previous one. The grid spacings were 
0.03125M/4M for the low resolution run, 0.02083M/2.66667M for the medium res- 
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olution run and 0.01563M/2M for the high resolution run. Here the numbers refer 
to the resolution on the hnest/coarsest grids respectively. The outer boundary of the 
domain is located at 200M and we impose periodic boundary conditions for simplicity. 
This puts an upper bound on the time up to which we can evolve the system before 
boundary effects influence physical observables. 

In Figure 8 (top) we display the real part of the i = 2, m = 0 mode of rT 4 extracted 
on a sphere of radius R = 60M using 4th order interpolation. We use 320 grid points 
in both the polar and azimuthal directions on the extraction sphere. Following [13], 
we test convergence by comparing a physical quantity T at different resolutions. The 
convergence is of order Q if for a set of grid spacings hi, ^ 2 , h^, the differences between 
the numerically computed physical quantity T at successive resolutions satisfy 

^ h^-h^ 

- 'i'hs hf - hf 

With the resolutions used in these runs, assuming 4**^ order convergence the above 
factor is ~ 5.953, whilst assuming 3'"'^ order convergence the factor is ~ 4.115. 

The gravitational wave content of the superimposed initial data is reflected in the 
non-zero initial signal. The collision of the two black holes takes place at f ~ 50 on this 
plot, so the signal before this collision time should be regarded as unphysical. As can 
be seen in this plot, the results for the two higher resolutions are indistinguishable on 
the scale employed here, whilst the lowest resolution shows a very slight drift towards 
later times, but is still in very good agreement. The bottom plot in Figure 8 shows 
the absolute value of the difference between computed low and medium resolution 
(solid blue), medium and high resolution (solid black), and this latter curve scaled 
up by the convergence factor assuming 3’"'^ (dotted orange) and 4**^ (dotted red) order 
convergence. This plot shows that in the highly dynamical stages of the evolution, when 
there is a lot of regridding and the boxes move around the domain, the convergence is 
closer to 3’"'^ order. On the other hand, when the system has nearly settled, and hence 
the boxes do not move much, the convergence order is closer to 4. We can explain this 
loss of convergence due to regridding because in the interpolation used in GRChombo 
only the values of the functions are matched across levels, but not their derivatives. 
We hope to improve this aspect of the code in the future. 

4.4 Choptuik scalar field collapse 

We now test the scalar held part of the code, by simulating the Choptuik scalar held 
collapse as described in [70] and illustrated in Figure 9. The referenced description 
is for a 1-1-1 simulation which is evolved using a constrained evolution, such that the 
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- Low resolution 

- Medium resolution 

- High resolution 
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Figure 8: Convergence test of head on collision. Top: The real part of the £ = 2, m = 0 
mode of r'I '4 on the sphere of radius R = 60M. Bottom: Differences in the real part of the 
i = 2, m = 0 mode of r'I ‘4 between three different resolutions. We also show the data rescaled 
by a factor consistent with either third (x4.11) or fourth (x5.64) order convergence. 


lapse a and the single degree of freedom for the metric, A, are both solved for on each 
slice using ODEs obtained from the constraint equations. The only degrees of freedom 
which are truly evolved are those of the held, 0, 4/ and If. 

Our evolution is carried out using the full 3 + 1 BSSN equations, without assuming 
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DB: KEx000000.3d.hdf6 dB: KEx000010.3d.hdf6 

Cycle: 0 Time:0 Cycle: 10 Tlme:5 



DB; Mlx000120,3d.hdf5 
Cycle: 120 Time;15 



Figure 9: In Choptuik scalar field collapse, the initial specially symmetric configuration 
in the first figure (which shows the values on a slice perpendicular to the z axis) collapses, 
splitting into an ingoing and an outgoing wave as seen in the second image. If the amplitude 
of the initial perturbation is greater than a certain critical value, the ingoing wave will result 
in the formation of a black hole, as seen from the output of the apparent horizon finder in 
the third figure, which shows that an apparent horizon with a mass of about 0.25 has formed 
by t = 15. 


or adapting coordinates to spherical symmetry. We are able to replicate the results 
obtained in [70], subject to some minor differences due to the fact that we evolve with 
the puncture gauge rather than according to the maximal slicing constraint equation, 
see Figures 23 and 24, which can be found in Appendix A. Videos of the results can 
be viewed via our website at http://grchombo.github.io. 

We see that GRChombo can accurately evolve the held prohle in the presence of 
gravity, and copes with the collapse of the supercritical case into a singularity, without 
code crash. For the subcritical cases we see that the held disperses as expected. 

4.5 MPI Scaling properties 

We now turn to the performance aspects of GRChombo. Here we perform a number 
of scaling tests to show that our code can exploit the parallelism ohered in modern 
supercomputers to a reasonable extent. Whilst Chombo does have the capability to 
partially utilise threads through hybrid OpenMP routines, we will limit our attention 
to pure MPI mode in these tests, as we have found that this gives signihcantly smaller 
run-to-run performance variations. 

Our strong scaling test is performed using a head on binary black hole system. 
We set up Brill-Lindquist initial data for two static black holes of mass 0.5M, with a 
separation of 6M. Our overall computational domain is a box of size 160M, and at the 
coarsest level, we hx the total number of grid points to 320 in each direction, giving a 
grid spacing of 0.5M. The centre of mass of the system is at the centre of the domain. 
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For the mesh rehnement, we hx the total number of levels to six. The simulation is 
allowed to run up to the time of 2M. The bulk of this test was performed on the 
SuperMike-II cluster at the Louisiana State University. Each compute node consists 
of two 2.6GHz 8 -core Sandy Bridge Xeon processors, connected via a InhniBand QDR 
fabric. We hx the computational load across all jobs and vary the number of core count 
from 16 to 2048. Our data in Figure 10 shows excellent strong scaling up to 200 cores 
on this cluster. We continue to see a reasonable speedup up to around 1000 cores for 
this particular problem. 
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Figure 10: Strong scaling behaviour of GRChombo on the SuperMike-II cluster at the 
Louisiana State University. The code achieves excellent strong scaling up to 200 cores, and a 
useful scaling up to around 1000 cores. 


Of course, in a production environment, it is often desirable to use additional cores 
to be able to run a larger simulation, rather than to speed up a problem of hxed size. 
In this scenario, weak scaling behaviour is of interest. We begin at 1024 cores with 
an identical setup to that in the strong scaling test. We then scale up the number of 
grid points at the coarsest level proportional to the increase in core count up to 10240, 
whilst adjusting the tagging threshold in order to maintain the shape and size of the 
rehned regions. We also adjusted the time step size (i.e. the Courant factor) so that 
each simulation would reach the target stop time in the same number of steps. We use 
the Mira Blue Gene/Q cluster at the Argonne National Laboratory for this due to the 
larger number of cores available. Figure 11 shows a less-than-perfect scaling behaviour 
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in this setup, with the main bottleneck appearing in the regridding and box generation 
stages. We are working together with the developers of Chombo to improve this aspect 
of the code performance. It is worth noting, however, that even in its current state the 
code still shows a useful level of scalability: the wallclock time increases by less than 
2 x over the lOx increase in core count. 
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Figure 11: Weak scaling behaviour of GRChombo on the Mira Blue Gene/Q cluster at the 
Argonne National Laboratory over a lOx increase in core count. 


4.6 Performance comparison 

Lastly, we demonstrate that GRChombo’s performance on standard 3 + 1 black hole 
problems is comparable to that of an existing numerical relativity code. 

Our comparison target is the Lean code [18,19], a 3 + 1 numerical relativity code 
designed to evolve four and higher dimensional vacuum spacetimes. Lean is based 
on the Cactus computational toolkit [74] and realises moving-box mesh rehnement 
via the Carpet package [75,76], both of which are part of the open-source Einstein 
Toolkit [13,22]. Initial data is constructed either analytically or numerically by em¬ 
ploying the TwoPunctures spectral solver [77]. In order to track apparent horizons. 
Lean makes use of AHFinderDirect [78,79]. 

The GRChombo setup is identical to that in the strong scaling test as detailed in 
Sec. 4.5. The Lean code is subject to the limitation of Carpet, where successive levels 
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may only occur a collection of nested-box hierarchies, whose sizes are typically related 
by a power of two. In this case, we hrst £x boxes of side lengths 160, 80, 40 and 20M 
at the centre of the domain, encompassing both black holes, then fix further boxes 
of side lengths 5 and 2.5M centred at each of the black holes. During the evolution. 
Lean has the capability to track the black holes and move or merge the hner boxes as 
appropriate, however the shape and size of the boxes remain unchanged. The GRChombo 
code is not subject to this box structure limitation, and therefore we simply tune the 
regridding threshold so that the size of the hnest level matches that of the Lean setup. 
We make no attempt to match the sizes of the intermediate levels as this would defeat 
the spirit of fully-flexible AMR. 
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Figure 12: Runtime and scaling comparison between GRChombo (orange) and Lean (blue). 
The leftmost data points show disproportionately large wallclock times as the machine be¬ 
comes memory-limited at this core count. 
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Our comparison tests were performed on the COSMOS VIII shared memory facility. 
Both codes were executed on the same SGI UVIOOO machine, utilising up to 60 Nehalem 
EX 2.67GIIz CPUs with 6 cores per CPU, giving up to 360 cores in total. In all of these 
runs, we pin one MPI rank to each core and disabled all checkpointing activity since 
we wish to exclude I/O bottleneck. We allowed the simulation to run up to coordinate 
time t = 2, and measured the wall-clock time taken to execute the time evolution 
portion of the code (i.e. we excluded the time spent during initial setup). 
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Within the range of 150-360 cores, both GRChombo and Lean exhibit similar perfor¬ 
mance and strong scaling characteristics (figure 12). Below 150 cores, we cannot mean¬ 
ingfully test the strong scaling behaviour as the machine becomes memory-limited. We 
have not performed this comparison on a larger cluster due to the lack of resource 
availability, but we have no reason to expect any significant difference provided that 
the problem size is also scaled up appropriately. Having said this, we believe that a 
framework like Cactus probably remains the better choice when it comes to these stan¬ 
dard problems, owing to the wealth of existing tools and resources and a more mature 
community of users. Instead, we intend for GRChombo to be complementary to existing 
numerical relativity codes in order to open up new avenues of research by enabling a 
wider range of problems to be tackled at a feasible level of resources (see section 3 for 
examples of such problems). 

5. Discussion 

In this paper, we introduced and described GRChombo, a new multi-purpose numerical 
relativity code built using the Chombo framework. It is a 3 -|- ID finite difference 
code based on the BSSN/CCZ4 evolution scheme. It supports Berger-Collela type 
AMR evolution with Berger-Rigoutsos block structured grid generation, and is fully 
parallelized via the Message Passing Interface, and time evolution is via standard 4-th 
order Runge-Kutta time-stepping. 

We illustrated some areas of physics that can potentially benefit from this code, 
such as multiple black hole mergers and scalar field collapse. Such fields require a code 
which adapts to changes in the range and location of scales at different points in space 
and time in the simulation. We emphasise that setting the initial conditions for these 
mergers are trivial - GRChombo automatically remeshes the grid given a set of analytic 
initial conditions without further user intervention. 

We showed that GRChombo successfully passes the standard “Apples with Apples” 
tests^. In addition to these tests, we evolved standard single black hole spacetimes 
(Schwarzschild and Kerr) and showed that it is stable to more than T = lOOOOM. 
Using the moving puncture gauge, we also show that GRChombo stably evolves the 
merger of two and three black holes in inspiral and head-on collisions. We simulated the 
supercritical collapse of a scalar field configuration, and found that it forms a black hole 
as expected, to show that the code supports non vacuum spacetimes. Finally we tested 

®We note that we perform as well as any other BSSN code in the Gowdy wave test as we expected 
in a pre-determined gauge. Although, [71] managed to achieve long-term evolution by considering 
different gauge conditions. 
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the MPI scaling properties of the code, both strongly and weakly, and compared this 
with an alternative numerical relativity code based on the popular Cactus framework. 

Nevertheless, despite its power, the AMR capability of GRChombo has to be treated 
with care. As we mentioned earlier, coarse-hne boundaries could be a signihcant source 
of inaccuracy, even though the Hamiltonian constraint may still be kept under control. 
A way to reduce coarse-hne boundary errors is to introduce conservative rehuxing 
during interlevel operations. Although rehuxing requires signihcant overhead, we intend 
to implement it in the next iteration of GRChombo. The litmus test for accuracy of 
GRChombo is its ability to make accurate predictions of outgoing gravitational wave¬ 
forms. We leave this, and the introduction of a set of “best practices” for the use of 
AMR in general relativistic systems, for a follow-up work. 
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A. Results from code tests 


In this appendix we collect the hgures of the code tests described in Sec. 4. 




(a) Hamiltonian constraint 


(b) Deviation in 


Figure 13: Robust stability test for both the BSSN and the CCZ4 codes, with resolutions 
p = 2,4 respectively. Left: time evolution of the Lqo norm of the Hamiltonian constraint. 
Right: deviation of jxx from 1. Neither norm grows at an increasing rate with increasing 
resolution, and so the test is passed. 
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(b) Error in Qyy at T = 1000 


Figure 14: Linear wave test. Left: analytical solution and the evolved gyy component of the 
metric at T = 1000 at resolution p = 4, but the two are indistinguishable. Right: absolute 
value of the error across the grid at T = 1000, from which we can see more easily that some 
small errors in the magnitude and phase have been introduced. 


42 






t [ct] 

Figure 15: Gauge wave test. The increase in the Loo norm of the Hamiltonian constraint 
means that the BSSN code only remains stable for less than 50 timesteps. Undamped (u) 
CCZ4, i.e. CCZ4 with ki = 0, performs similarly. Damped (d) CCZ4 with ki = 1 is stable 
for the full 1000 timesteps. The test was performed with initial amplitude of A = 0.1, Kreiss- 
Oliger dissipation coeffecient of u = 0.1 and a resolution of p = 4. 
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(a) BSSN Lapse a (collapsing) (b) Hamiltonian constraint (collapsing) 

Figure 16: Gowdy wave test, collapsing. Left-, minimum value of the lapse a across the grid 
as the spacetime is evolved towards the singularity. As expected, the harmonic gauge causes 
the evolution to “slow down” as the singularity is approached. Right-, evolution of the Lqo 
norm of the Hamiltonian constraint for two resolutions. The test reaches T = 1000 crossing 
times without crashing. 



(a) K (expanding) 



(b) Hamiltonian constraint (expanding) 


Figure 17: Gowdy wave test, expanding. Left: trace of the extrinsic curvature K as the 
Gowdy wave spacetime is evolved in the collapsing direction. This correctly asymptotes to 
zero as the spacetime expands, but becomes unstable at around t = 30. Right: evolution of 
the Loo norm of the Hamiltonian constraint for two different resolutions. 
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(a) Convergence (expanding) 


(b) Convergence (collapsing) 


Figure 18: Gowdy wave test, convergence. The ratio of the L^o norm of the Hamiltonian 
constraint for the resolutions p = 4 and p = 2 is shown, for the expanding and collapsing 
directions for the BSSN and CCZ4 codes. A value of 16 indicates 4th order convergence, 
which is demonstrated by the codes initially, although lost at later times by BSSN. 



(a) Hamiltonian constraint 
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(b) ADM quantities 


ADM Mass 
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Figure 19: Schwarzschild black hole simulations. Left: Evolution of the norm of the 
Hamiltonian constraint up to t = lOOOOM, showing long term stability. Right: ADM Mass, 
angular momentum and linear momentum (in the x direction) during the initial stages of the 
evolution. These quantities remain approximately constant. 
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(a) Hamiltonian constraint 
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(b) ADM quantities 
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Figure 20: Kerr black hole simulations. Left: Evolution of the norm of the Hamiltonian 
constraint. Right: Components of the angular momentum and mass of the Kerr black hole 
during the evolution. The ADM quantities remain constant. 
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(b) ADM linear momentum 


Figure 21: Boosted black hole simulations. Left: Evolution of the L^ norm of the Hamil¬ 
tonian constraint. Right: Components of the ADM linear momentum during the evolution. 
They remain constant. 
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Figure 22: L?' norm of the Hamiltonian constraint across the whole domain for the binary 
black merger. The constraint violation remains bounded throughout the evolution, which 
includes the merger and ring-down phases. 




(a) Subcritical profiles of </) 


(b) Supercritical profiles of cj) 


Figure 23: Choptuik scalar field collapse. The profiles shown for the fields at t = 0, 5 and 
20 differ from those in [70] due to the different gauge conditions used. In the supercritical 
case we show the snapshot at t = 10 rather than 20 as this is the point at which the evolution 
is frozen in the gauge choice in [70]. In the puncture gauge the evolution of the region within 
the event horizon continues and the result is that the large spike in the held effectively falls 
into the puncture, resulting in a zero held value at the centre of the coordinate grid. 
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(a) Subcritical lapse (b) Supercritical lapse 

Figure 24: Choptuik scalar field collapse. The values of the lapse at the centre of the grid 
are given. It can be seen the the profiles are very similar to those obtained by Alcubierre 
in [70], and that the one for the supercritical case shows the characteristic collapse of the 
lapse which is symptomatic of black hole formation. 
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